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Abstract 

We consider a single particle spectrum as given by the eigenvalues of the 
Wigner-Dyson ensembles of random matrices, and fill consecutive single par- 
ticle levels with n fermions. Assuming that the fermions are non-interacting, 
we show that the distribution of the total energy is Gaussian and its variance 
grows as n^logn in the large-n limit. Next to leading order corrections are 
also computed. Some related quantities are discussed, in particular the near- 
est neighbor spacing autocorrelation function. Canonical and grand canonical 
approaches are considered and compared in detail. A semiclassical formula 
describing, as a function of n, a non-universal behavior of the variance of the 
total energy starting at a critical number of particles is also obtained. It is 
illustrated with the particular case of single particle energies given by the 
imaginary part of the zeros of the Riemann zeta function on the critical line. 
Submitted for publication to Physica D - April 27, 1998. 
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I. INTRODUCTION 



Consider a system of n interacting fermions. The ground state total energy may often be 
well approximated by macroscopic methods, a prototype of these being the Thomas-Fermi 
approximation. Aside from the well-known expansion containing the volume energy, surface 
energy, etc, these methods can be further improved by the inclusion of additional systematic 
quantum effects, like for example shell effects related to some symmetries of the system (see 
for instance [Q). Once all these contributions have been taken into account, one would like to 
have some estimate of the remaining "irreducible" discrepancies between the calculated and 
measured total energies due to non-systematic effects. These discrepancies may be viewed 
as fluctuations observable as some parameters of the system are varied, like for example the 
total number of particles or any other relevant external parameter. 

On the other hand, it has been established that statistical properties of systems whose 
classical analog is chaotic are well described by eigenvalues of ensembles of random matrices 
(see for instance 0). It seems therefore reasonable to suggest that an extreme estimate of 
the above mentioned fluctuations may consist of modeling the single particle spectrum by the 
spectrum of a random matrix of the Wigner-Dyson type. This may be a crude approximation 
because the statistical properties of the low-lying energy levels of the single particle sequence 
may not be well described by a random matrix approach. Another shortcoming of such a 
model is connected to the fact that, for simplicity, we fix to a constant the average density 
of single particle states, and ignore its variations with energy. In its simplest version, the 
model may apply either when only properties around the Fermi level are considered or when 
particles are enclosed in a 2-dimensional box. However, inclusion of a variation (increase) of 
the density of single particle energies can be taken care of. This model may also be relevant 
in studying the energy distribution of systems of interacting fermions trapped in a chaotic 
enclosure, like for example the variations of the total energy of n electrons contained in a 
quantum dot as the shape of the dot is varied by some external potential or the variation 
produced when varying a magnetic field. 

To be specific, consider an ordered sequence of energy levels xq < xi < . . . < Xi < . . .. 
Let Si, S2, S3, ... be the distances between consecutive eigenvalues. 



We assume that the sequence is stationary and we set the mean spacing D to one. We 
are interested in the statistical properties of the "ground state" energy of a system of 
n non-interacting fermions resulting from adding the "single particle" energies Xi of the 
i = 1, 2, . . . , n occupied states 
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where xq has been taken as the energy origin. Consider now an ensemble of single particle 
spectra which for the moment we let unspecified. We expect under very general conditions 
the probability distribution of Xf to be Gaussian in the large n limit. This, by the central 
limit theorem, is obviously true for a set of uncorrelated levels. For (correlated) eigenvalues 
of Gaussian ensembles of random matrices this has also shown to be true, in fact in a more 
general context . As a numerical illustration. Fig. |I] displays the probability distribution 
of Xt for n = 300 computed for a Gaussian orthogonal ensemble of random matrices. Our 
purpose is to derive expressions for the average and the variance of Xt, thereby specifying 
completely the "ground-state" energy distribution. 

The first physical situation we consider is when the number n oi particles is fixed. This 
case, which we call ^^canonicaT by analogy to conventional statistical mechanics, is treated in 
section II. From the mathematical point of view it is related to the autocorrelation function 
of consecutive spacings, a poorly known quantity which is not a two-point measure. Another 
possibility is to consider the typical fluctuations of the total energy of the occupied levels 
contained in an energy interval of given length. In this case the number n of levels fluctuates 
from sample to sample around its mean value. This alternative ^^grand canonical^ problem 
will be treated and compared to the canonical one in section III. The grand canonical case 
has the advantage of being expressible in terms of the spectral two-point correlation function. 
As expected we find that both approaches give the same asymptotic increase of the variance 
of the total energy (proportional to ra^logn in the large-n limit). 

The results of sections II and III correspond to the "universal regime" as given by ran- 
dom matrix theory. In section IV we use standard semiclassical techniques to show that 
generically for ballistic - as opposed to diffusive - systems there is a saturation effect, 
namely a crossover from the log n to a growth of the variance of the total energy as a 
function of n. This new regime holds when the number of particles is larger than a system- 
dependent critical value and is accompanied by non-universal oscillations around the mean 
growth for which we give an explicit expression. To illustrate this point with a particular 
example, we have considered the unphysical but interesting in its own and explicitly com- 
putable case where the single particle energies are given by the imaginary part of the zeros 
of the Riemann zeta function on the critical line. Section V summarizes the results and 
gives some perspectives. 



Statistical averages over the ensemble are denoted by the symbol (.). The mean value 



II. CANONICAL VARIANCE 



{Xt) is 




n n 




(2.1) 



i=l i=l 



2 



where the mean level spacing has been set equal to one. 

Let us now consider the variance of Xt defined, from Eg . (|1 . 1|) 



by 



n 



/S?{n) ^ {X't) - {Xt)' = ^ (n - 2 + l)(n - J + 1) (s, s,) 
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ran 

i=i i^j=i 

The main quantity of interest in this expression is the autocorrelation function I{k) of 
spacings of consecutive levels, 

I{k) = {siSi+k) /c = l,..,n-l. 

The stationarity of the spectrum implies that I{k) depends only on the relative index k and 
that {sf) = (s^), for all i. Performing the sums over the index i 

, n, , n(n + l)(2n + 1) , 1 9 , ^9 1 ^ ^, , ^, , , , 

A\n) = ^ ^-{s^) - -n' (n + iy + -J2 Hn^j) /(j), n>2, (2.2) 

where 

F{n,j) = f - (3^2 + 3n+ l)j + n {n + 1) (2n + 1) . 

For n = 1 the variance takes the value A^(l) = (s^) — 1. All the non-trivial information is 
contained in the autocorrelation of spacings; a detailed study of this function will be made 
in the next subsection. The extreme case of a uniform uncorrelated sequence of levels is 
particularly simple since I{j) = 1 for all j, and it follows from Eq.( |2.2| ) that 

A^-((.^)„-l)(^4 + 9). (2.3) 

When the nearest neighbor spacing distribution of the uncorrelated sequence is -P(s) = 
exp(— s) then in Eq.( p.3| ) is equal to 2. 

A. The autocorrelation function of spacings 

For Gaussian ensembles (GE) of random matrices, the quantities entering Eq.( p.2|) can 
be expressed as 

roo 

(s^) = 2 / ^(0,s) ds , 
Jo 

I{j)= / E{j,s)ds, j>l, (2.4) 





where E{j, s) is the probability that a randomly chosen interval of length s contains exactly 
j eigenvalues |0. The functions E{j,s) can be written in terms of an infinite product of 
the eigenvalues of an integral equation involving spheroidal functions and only numerical 
tables of these functions for low values of j and a limited range of s exist. Therefore the 
integrals in ( |2.4|) cannot be computed analytically and explicit expressions for the I{j) are 
not available. Some numerical estimates of the integrals (273), to be used in what follows, 
are g 

(s^) = { 1.180, P = 2 (2.5) 
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where P =1,2 and 4 denote the three types of GE, orthogonal (GOE), unitary (GUE) and 
symplectic (GSE) respectively. Numerical estimates of exist for low values of j. For 
instance, for j = 1 H 



1(1) 



0.922, 
0.944, 
0.964, 



P = l 
P = 2 
/3 = 4. 



(2.6) 



For our purpose the functions I{j) for arbitrary j are needed and the following ansatz 

Hj) = 1 - v/f - A/j' - a/f (2.7) 

will be made. The constants in Eg. ( ^.7[ ) are, in principle, dependent on the symmetry class 
of the GE considered {j3 = 1, 2, 4). Eq. (|2.7|) is one of the basic equations of this paper. The 
simple functional dependence on j is consistent with the requirement that two sufficiently 
far apart consecutive spacings should be uncorrelated. Moreover, general considerations 
suggest that odd powers of j may be excluded. As we will now show, the three parameters 
entering Eq. ( p.7| ) may be uniquely determined by requiring the correct asymptotic behavior 
of the fluctuations in the length of an interval containing a fixed number of levels. 

Consider for that purpose the statistical fiuctuations in the length 5* of a spacing made 

n 

of n consecutive nearest-neighbor spacings, 5* = ^ . The spacing variance o"^ of S can be 

i=l 

expressed in terms of the I{j) in the following manner 



n-l 



(n)^ {S')-{Sf=n{{s')-n)+2Y.in-j)I{j), n>2 



(2.8) 



while cr^(l) = (s^) — 1. We use a non-conventional but more natural notation, namely 
(T^(n) denotes the variance of the distribution of the sum of n consecutive nearest neighbor 
spacings (n = 1, 2, . . .). The conventional notation being (T^(n — 1), n = 1, 2, . . .. 
Replacing the ansatz (|2.7|) in (|2.^ ) we get 

a^{n)= ((s^)-l) n + 2ri [^(n) + 7] + 2 r/ n [^(^^(n) - C(2)] + A [¥^\n) + 2 ((3,) 



+ 



n A/3) \^^^\n) - 6 C(4)] + (a/12) \^^^\n) + 24 C(5) 



+ (n a/60) 



n) 



120 C(6) 



(2.9) 



n— 1 



where 7 is the Euler constant, \E'(n) = —7 + - is the DiGamma function and ^^"''\z) = 

^"^(z) are the PolyGamma functions 0. The large-n behavior of (T^(n) follows from the 
asymptotics of the vl/^^^'s 0. Only \1/ and "^^^^ will be needed, the other not contributing 
at the order we are working 



vl>(n) = logn - ^ - ^ + 0(l/n^) 



(2.10) 
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Then 



= An + 2r] logn + B+ 0{l/n^) , (2.11) 

with 

A^(,2)_i_2ryC(2) -2AC(4)-2aC(6), 
5 = 2(7 + l)r/ + 2AC(3) + 2aC(5). 

Would one have retained a 1/j term in the ansatz, the leading order term would have been 
of order n log n. 

A different quantity, closely related to the spacing variance ^^(n), is the number variance 
T?{L) defined as the variance of the number of levels contained in an interval of given length 
L taken at random. This quantity is known analytically for the three Gaussian ensembles 
and has in the large-L limit the asymptotic behavior 0] 

S2(L) = -^logL + C + 0(l/L) , (2.12) 

where 

C = 2/(/?7r2){7 + log(27r) + 1 + /?(/? - 2)7rV8 + [log2 - l-n^l^^d^^^] . (2.13) 

Consider now more carefully the relation between the spacing and the number variances. 
The values that the number variance takes at integer values of its argument are, for large n, 
related to the spacing variance through |§| 

T?{L = n)-a\n) = llQ . (2.14) 

This equation indicates that: (i) as expected from their definition, the leading order terms of 
both quantities coincide and, (ii) the next-to-leading order terms differ by a constant (1/6), 
irrespective of the symmetry class j3. 

Eq.( |2.14j ) completely determines all the unknown coefficients in Eq.( p^.7D . In fact, con- 
sistency with Eq. ( p.l4| ) imposes the following conditions in Eq. (p.ll| ) 



1 



A = 0, (2.15) 
C - = 1/6 . 

The remaining free parameters A and a are adjusted to satisfy the second and third condi- 
tions 

^ (s^) - 1 - 2 C(2) _ 06) 

2C(4) C(4) ' 

- 2[C(3)C(6W(4)C(5)] {^-^(^^^^-^n^ i[(^Vl-2.C(2)]) . (2.16) 
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Using Eq.(g]5]) these two constants take the values A = 0.02027, 0.03013, 0.00997 and 
a = -0.04483, -0.02550, 4 x 10"^ for /? = 1,2 and 4, respectively. 

In order to test the accuracy of Eq. ( p.7| ) we have computed numerically the autocorrela- 
tion function for orthogonal as well as unitary random matrices. As can be seen from 
Fig. the agreement between the ansatz and the numerical simulations for (3 = 2 is very 
good. Because the overall agreement is of similar quality we do not show the analogous 
curve for GOE. As an additional quantitative test let us mention that when expression ( p.7| ) 
is evaluated at j = 1 it reproduces the values in Eq. (|2.6|) with a 1 x 10~^ error for any f5. 
This is remarkable if one keeps in mind that only asymptotic information has been used 
when determining the parameters in the ansatz ( |2.7| ). 

In a previous study of the autocorrelation function of the spacings Odlyzko P has quoted 
a conjecture proposed by Dyson: J(j) 1 — rj/j"^ (see also chapter 16 in 0). This ansatz 
coincides to lowest order with Eq. ( p.7|) . Aside from considerations related to accuracy, let us 
emphasize that the inclusion of a term of higher order (like 1/j^ or in ( |2.7| ) is essential 
in order to obtain the correct asymptotic behavior for a'^{n) as well as, as we shall see later, 
for A^(n). This is so because these higher order terms ensure the vanishing of the linear 
term in Eq.( |2.1l| ). 



Eq. ( |2.7| ) is consistent with known relationships satisfied by the autocorrelation function. 
For example, the sum rule 

oo 

2E[Aj)-1]+^'(1) = 
i=i 

valid for GE of random matrices []TU[ is equivalent to the condition A = in Eq.( |2.15| ) (i.e., 
the vanishing of the linear term of a'^{n) |11]). 



B. Variance of the total energy 



We have now the necessary ingredients to compute the variance of the total energy. 
Replacing in Eq.( p.2|) the ansatz (|2.7| ) and using the definition of the DiGamma and 
PolyGamma functions it follows that 
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n 



2n^ 



1 
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for n > 2. Using again the asymptotic expressions (|2.10|) for \I'(n) and \I'*^^)(n), as well as 
those of ^'(^^(n) and ^'(^^(n) we finally get 



A (n) = 7] n log n+-[C — 7] — - ]n + r] n log ''T' + 7: [C — - ) n 



1 



1 



6 



7 



1 / ^^, 1 

6 b 



A-2aC(3) 



+ 0{l/n). (2.17) 



Notice that, as it happened for a'^{n) with the terms of order nlogn and n, the functional 



form of the ansatz ( p.7|) together with the method by which the parameters r], a and A have 
been determined insure the vanishing of the terms of order log n and of the variance of 
the total energy. That the remaining leading order behavior ra^ logn in Eg. ( p. 171) is correct is 
further confirmed by the fact that it coincides with the asymptotic result obtained from the 
linear statistic theory (see section III). This behavior is in contrast with the faster growth 
of order of an uncorrelated spectrum (cf Eq.(2J)). The difference is of course due to the 
rigid nature of the GE spectra. 



Fig. |] shows for GOE the comparison of numerical results with Eq. (|2.17|) for the variance 
of the total energy as a function of n. The overall agreement is good (for n 50 the relative 
error is less than 2%). The theoretical curve is sensitive to the precise value of (s^) which 
enters in the definition of A and a, and is known only up to the third digit (cf Eq.( p.5D ). For 
comparison we have also plotted in Fig. ^ the leading-order term of A^, which clearly fails 
to reproduce accurately the numerical results in the interval of n displayed (the error is of 
the order of 10% at n ~ 50). 

In the previous analysis we have taken a single particle level (xq in (|1 . 1| ) ) as the reference 
energy. In some physical applications it may be more appropriate to measure energies with 
respect to an arbitrary origin, which will not coincide with a single particle level but rather 
will lie in a random position between two of them. As shown in Appendix A when this 
construction is adopted, namely filling n successive levels located just above this origin, 
Eq.(p.l|) giving the mean value of the energy is modified as follows n(n + l)/2 — >• n^/2, 
whereas for the variance of the energy one has 



n , 



n] 



(2.18) 



Notice that both sides of (|2.1^ 
for higher order terms. 



have the same leading order, though different coefficients 



III. GRAND CANONICAL VARIANCE 

Here we treat a slightly different problem with respect to the previous section. Instead of 
fixing the number of occupied levels, we consider an energy interval of given length 2L and 
compute the energy variance for the single particle energies xi, . . . ,x„ contained in it (to 
follow conventional notation, in this section we will denote the length of the energy interval 
by 2L, whereas in the previous section we were denoting it by L). The length 2L of the 
energy interval is now kept fixed but the number n of levels contained in it may now vary 
from sample to sample. This is a particular example of a more general class of problems 
usually called "linear statistic" in random matrix theory [0, dealing with the distribution 
of a variable W of the form 

n 

W = Y.f{x,), (3.1) 

i=l 
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where / is an arbitrary function. The variance of W was given in 



12' 



Vw= [ t"" dx dy [6{x -y)- Y^i^x, y)] f{xD)f{yD) (3.2) 

or, ahernatively, in the Fourier space 

1 r°° 

Vw = - dk[l-h{Dk)]^{k)^{-k) . (3.3) 

U J— oo 

Here D is the mean level spacing (assumed to be constant; it will be set to one at the end), 
Y2 the two-level cluster function, h{k) = J°°^Y2{x) exp[2i'n'kx] dx the form factor and ip{k) 
the Fourier transform of f{x) 

ip{k) = J f(x) exp[—2iTckx] dx 

(/(x) is assumed to be zero outside the interval [—L,L]). 

To compute the number variance via Eg. ( |3.3| ), one must choose f{x) = 1, as Dyson 
and Mehta did. For our purpose, namely to compute the variance of the total energy, f{x) 
should be chosen as follows 

\ X + L -L<x<L 
= I \x\>L ^3.4) 

whose transform is 

^(k) = —^[2nLk cos(27iLk) - sin(27rLA;)] + — sin(27rLA;) . (3.5) 
2'K^k^ nk 

In Eq. (|3.4|) we have chosen to measure the energy with respect to the lower (or equivalently 
upper) end of the energy interval considered. This, which is somewhat natural and analogous 
to Eq. ( |1.1|) , is by no means compulsory and the reference may be taken at any given point 
of the energy axis. The variance of the total energy depends on the location of the reference 
point, the minimum being obtained when it is at the center of the interval. A general relation 
connecting the variance at different reference points is given in Appendix B. 

Let us now consider the case of a large (average) number n of levels contained in the 
interval. Though we are using for notational simplicity the same symbol as in the canon- 
ical case, for grand canonical expressions n will mean the expectation value (n) = 2L/D. 
Therefore, in contrast to the canonical case, it takes continuous values. Following Ref. ||12|| , 
we split the integral in Eq. (p. 3D in two parts 



2 /■£ 2 
Ve = — ip{k)ip{-k) [1 - b{Dk)] dk + — ip{k)ip{-k) [1 - h{Dk)] dk (3.6) 
D Jo D Je 

where we have used the parity of the integrand. The parameter e is chosen such that 
l/L-Ce<Cl/Z). For k < e, Dk is always much smaller than one and the form factor can 
be approximated by P 

biDk) ^l-'^\k\ 
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and the first term in the r.h.s. of Eg. ( p.6|) (denoted V^""^) can be written 

dx , 



1 2 cos X sin x sin^ x 

+ 



X x^ 



^ /37r2 Jo 

which when integrated and taken in the hmit Le ^ 1 gives 

Ve^ ^ [log(27rLe) + 7 + log 2 - 1/2] . (3.7) 

On the other hand, for k > ewe have kL ^ 1 and the function ip{k) shows rapid oscillations 
as compared to the variations of the form factor. We then replace the product ip{k)ip{-~k) 
in the second term of the r.h.s. of Eq.( p.6|) (denoted V^'') by its average value over k 



^ Tl^ JDe X2 ATT^Joe 

The second integral in the latter equation is of order l/(eL)2 with respect to the first one 
and therefore can be neglected. Evaluating the first integral we get 

—47' 2 7-2 
^ ^ [iog{De) - 1] + (/5 - 2)- . (3.8) 

Collecting the two terms we obtain the final expression for the variance which we express 
in terms of the average number n of levels contained in the interval and the constant C (cf 
( p.l3| )) (we put moreover D = 1) 

VE{n) ^ T] n'^\ogn+^{C - r])n^ . (3.9) 

Let us now compare grand canonical and canonical results. When the sequence consid- 
ered is large (n 3> 1), we expect the energy variance computed for a fixed number n of 
occupied levels (A^, canonical) to be equal to the energy variance when the fixed length of 
the energy interval contains on the average n levels {Ve, grand canonical). The comparison 
of Eqs. (|3.9| ) and ( p. 171) shows that indeed to leading order Ve = (the same leading order 



behavior has been obtained in Ref. [0 for the particular case of GOE employing Dyson 
Brownian motion). However the correction terms are different. It is interesting to note 
that the next to leading order corrections differ, as in Eq.( p.l4D for the number and spacing 
variances, by a constant independent of (3 



VeM - A^(n) 1 



which follows from Eqs.( |2.l7D and ([ 

To conclude this subsection let us point out that exact expressions for the energy variance 
can be computed in the grand canonical approach for any value of n without assuming that 
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it is large. For example for the unitary case (/3 
Eg .( 13 .51) and the fact that 



b{k) 



1- 




k 



2) from the definition Eg. ( |3.3| ), using 

(3.11) 



|< 1 
|> 1 



it follows that 



VEin) 



n 



2 r 



27r2 



47m, 



log 72 — Ci(27m) Si(27rn 



sin (27m) 



67m 



+ 



^ - ^) -«<27m) 



1 



127r%2 



2Tx'^n , . s 1 

+ 7 + log(27r) + - 



(3.12) 



where Ci(a;) and Si(x) are the Coslntegral and Sinlntegral functions 0, respectively. Ve(?t.) 
grows like for n — ^ (remember that here n is a continuous variable) and then switches 
to the \ogn behavior for larger values of n. 



IV. SEMICLASSICAL TREATMENT: NON-UNIVERSALITIES 

A. Chaotic Systems 

The strength and vitality of random matrix theories rely heavily on the universality of 
some of its predictions. As a counterpart it is poorly adapted to capture some system- 
specific features. However special tools have been developed for that purpose to study 
systems whose classical analog is chaotic. Indeed, one of the important achievements of 
semiclassical theories has been to determine the limits of validity of the universal regime as 



described by random matrix theories by including system-dependent corrections flj]. Our 
purpose now is to adapt these methods, which have been developed for the form factor and 
are thus applicable to the grand canonical case, to the study of the variance of the total 
energy. 

Our starting point will be Eg. (p.3|) written in a slightly different form 



2 r'-^ 

VW = dTK{r)^{r/D)^{-r/D) . (4.1) 

Here K{t) = 1 — 6(r) and we have used the parity of the integrand. For simplicity we 
restrict in the following to systems having no time reversal symmetry {(3 = 2). The results 
can be easily generalized to the other symmetry classes. 

In the previous section we computed the universal behavior of the variance of the total 
energy by inserting in ( [4.1| ) the random matrix spectral form factor (cf Eg. ( ^.llD for the 



unitary ensemble). It is possible however to give a more accurate description of the form 
factor based on semiclassical approximations. For systems having a classical analog the key 
ingredient is a formula expressing the spectral density p{E) = '^6{E — Ej) as a sum over 

j 

the classical periodic orbits 

-1 00 

piE) ^ {p{E))e + ^ E E Ap,r COS {rSjn) . (4.2) 

p r=l 
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Here {p{E))e = 1/D is the spectral average density of states, Ap^^ = ^p/y I det(M^ ~ 1)1 is 
an amplitude that depends on the period Tp and the stability matrix Mp of the primitive 
orbit p, r are the repetitions and Sp is the action We restrict our considerations to 

chaotic systems for which periodic orbits are isolated and unstable. 

Using Eg. (|4.2|) for the density of states Berry proposed the following form for K{t) 
for chaotic systems 



{D/h)J2Al5iT-rTp 



p,r 



T 
1 



r < r* 

T* <T <1 
r> 1 , 



(4.3) 



where h is Planck's constant and r = TD/h is the (rescaled) time measured in units of the 
Heisenberg time. The main physical ingredient entering this formula is a classical sum rule 



due to Hannay and Ozorio de Almeida which takes into account the exponential prolif- 
eration of long periodic orbits as well as their ergodicity. As a consequence of this as well as 
other semiclassical considerations the form factor for times r larger than a certain critical 
time T* becomes "universal" (i.e., coincides with random matrix theory). The situation 
is different for short periodic orbits which do not display this universality and whose con- 
tributions, explicitly written down in Eq.( ^.3| ), produce system-dependent deviations from 
random matrix theory for r < r*. Usually the parameter r* is chosen in order to satisfy 
■TVnin <^ T* <C 1, where Tmin = T^inD / h is the rescaled period of the shortest periodic orbit. 

The contribution of short periodic orbits to the variance of the total energy is obtained 
by replacing K{t) for r < r* in Eq. (|4.1|) , with the result 



^Ap^r ^{rTp/D) ip{-rTp/D) . 



(4.4) 



p,r 



Here Ap^r = h 'Tp / {D\J\ det(Mp — 1)|) is the usual amplitude Ap^r written in terms of the 
rescaled time Tp. 

For times r > r* we can repeat the same steps as in the previous section (the form factor 
coincides with GUE), with the important difference however that now there is a cutoff in 
the integral ( [4.1D at r = r*. The variance of the total energy may now be written in the 
following form 



(2), 



(4.5) 



where V^^^ is the GUE result (|31^) and 



27r2 



cos(27mT*) sin(27rrir*, ^ , 
log nr* + ^ + ^ - Ci 27mr* 

+ log(27r)+7-i 



A2 
^ p,r 



rTp<T* p 



COS vmrr, 



2TTnT* 



sm 7rnrT„ 



rmrTr, 



1 



{27mT*y 



+ sin [TcnrT„ 



(4.6) 
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(2) 

is the contribution to the variance due to the lower hmit at r = r* in the integral 
3]). Aside from the system-dependent parameter r*, its functional form is general. V^^^ 



is obtained from ( ^.41 ) using the explicit form ( |3.5| ) of the Fourier transform Lf and contains 
detailed system-specific information. 

The energy variance as a function of the number of particles exhibits now two different 
regimes. For n < 1/r*, and V^^^ are of order (nr*)^ and can therefore be neglected. 
The energy variance is then given by the random matrix expression. For n ~ 1/r* these two 
terms can no longer be neglected and there is a crossover to a different regime, not described 
by random matrix theory. This is clearly seen for n ^ 1/r* when there is an almost perfect 
cancellation between V^^^ and leading to 

7~)2 2 

VE{n)^—^\l-\ogT* + 0{l/nT*f]+vi\n) n > 1/r* . (4.7) 



From Eq.(4^) for V"^^'' we then see that now the variance instead of growing as logn like 
in random matrix theory, saturates and increases as k n^. The prefactor k is not constant 
but shows non-universal oscillations as a function of n described by the short-periodic-orbit 
contributions. In order to determine its average value over n we replace each oscillating 
term between curly brackets in Eq.( [4.6| ) by its mean value (which is equal to one). Then, 
denoting this average by an upper bar 

P D^n^^ 1 



^(3) _ ^ II' sr^ ^p,r ^ II sr^ 

~ 27r2/i2 r2r2 ~ 2ti^ ^ det(M; 



where the last sum extends over all orbits satisfying rTp < hr* /D. An estimate of V^'^ 
may be obtained by replacing this sum by an integral weighted by the density exp{hKsT)/T 
over the period T of the periodic orbits {hxs is the Kolmogorov- Sinai entropy) |T^. This 
approximation, which is not well justified in this regime of "short periodic orbits", will 
be shown however to give good results for the zeros of the Riemann zeta function (see 
next subsection). Then using the Hannay-Ozorio de Almeida sum rule | det(Mp — 1)| — >■ 
exp{hKsT) one has 



Inserting ( [4.8| ) in ( |4.7] ) we finally obtain 



VE{n) = —;- l-logr^i„ + 0(l/nr*)2 . (4.9) 

Thus in this approximation the average of the energy variance in the non-universal regime is 
determined by a single parameter, namely the rescaled period of the shortest periodic orbit. 
Not only is the average of the energy variance determined by the shortest period. For 



n ^ l/rmin the contribution of the term sm{TTnrTp) / (TrnrTp) in Eq.(4^) is negligible for any 
Tp. Then, in this limit, each oscillating term within curly brackets is again equal to one and 
therefore the energy variance itself is given by the r.h.s. of Eq. ( |4.9| ). In summary, in the 
non-universal regime the energy variance shows, superimposed to an n"^ growth, oscillations 
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whose amplitude is damped with increasing n. As the number of contributing terms in ( |4.6|) 
diminishes, the structure of the oscillations becomes more regular and eventually only the 
lower frequencies survive just before the complete extinction of the oscillations. 

It is worth noticing also that the damping of the oscillations is a remarkable peculiarity 
occurring only when the reference point used to measure the single particle energies is at 
one border of the energy interval considered. For example, setting the origin at the center of 
the interval, no damping is found and the amplitude of the fluctuations remains constant as 
a function of n. This latter behavior with non-universal oscillations of constant amplitude is 
similar to what happens for the other two-point measures mainly considered so far, namely 
the number variance and the Dyson-Mehta least square statistic [p!^p!7[ . 

It would be valuable to test the previous results (and in particular Eg. ( [4. 51) ) as a function 
of n for a real physical system, like for example by locating of the order of 100 particles 
inside a billiard of nanometric size as in the coulomb blockade experiments in quantum dots 
T^ . Another possibility, not related to experiments but which can be explicitly computed. 



is the Riemann zeta function. 



B. Application to the Riemann zeta function 

This function constitutes a paradigm in the field of quantum chaos. There are many 
reasons for that. On the one hand the density of zeros of C(l/2 + iE) as a function of E 
may be expressed through the equation 

pdE) = {pdE))E-lT.f:y^^os{rElogp) . (4.10) 

The sum here goes over all the prime numbers p = 2,3,5,... and {P(^{E))e = 1/D = 
(l/27r) log{E/2'ir) (we are systematically ignoring here problems related to the convergence 
of the series). Eq.( |4.10| ) has the same structure as Eq.(^]^) and this suggests a dynamical 
interpretation of it. On the other hand, it has been shown numerically |^ and proved in some 
cases and demonstrated by heuristic arguments in others [p!^,p|, p!^pO| that many statistical 
properties of the critical zeros of that function coincide with the GUE results of random 
matrix theory. 

Because of the resemblance of the properties of the critical zeros of the Riemann zeta 
function with those of the eigenvalues of a real chaotic system, and because of the possibility 
of performing explicit computations for that function, our purpose now is to use the imagi- 
nary part of the zeros as the "single particle spectrum" of some unknown physical system. 
As a first test of the results obtained in this paper we have considered the nearest neighbor 
spacing autocorrelation function I{j) for a set of 50 000 zeros located around the lO^^th zero 
of C(l/2 -l-iE) (the set starts at Eo = 267653395648.8475 . . .) computed by Odlyzko 0. The 
values of for low j, shown in Fig. H, are in good agreement with those obtained from 
the ansatz ( p.7| ). Furthermore, we have computed for the same set of zeros the variance of 
the total energy as a function of the number of occupied levels. Before presenting the results 



we need to compute explicitly the third term in the r.h.s. of Eq.([4.5|), denoted V^l. An 



analogous contribution for the number variance of the Riemann zeta function was already 



considered in 
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By comparing Eqs.( |4.2| ) and ([4.10|) one can see by analogy that the period of the pe- 
riodic orbits should be identified to the logarithm of the prime numbers, Tp = logp, and 
moreover h = 1. Then the rescaled period is Tp = logp/ log(i?o/2vr) and, from ( [4. 101 ), 
-^p,r = ~27rrp/(Z}p''"/^). Then, as for any dynamical system, it is expected that K{t) be- 
haves in a universal manner for t > t*, while "small" prime numbers should contribute with 



non- universal terms for r < r* . From the expression of 
mentioned identifications we get 



(3) 



in Eq. ( |4.6D and using the above 



V, 



(3) 



rTp<T 



27r2 



cos[7rnrTp) — 



sm[7inrTp 
nnrTn 



+ sin^ (TmrTp) 



(4.11) 



For the zeros we are considering we have Tmin = log2/log(£'o/27r) ^ 0.028. We then choose 
r* = 0.2 for the numerical comparisons, which leads to a maximum prime number Pmax = 131 
in the sum in Eq.( ^.ll| ). 

In order to magnify the effect of non-universal corrections we have plotted in Fig. ^ 
the normalized quantity VE{n)/{n'^ /27t'^) as a function of n (we set D = 1). The solid line 
represents the theoretical prediction ( |4.5| ) with the third term of the r.h.s. given by Eq. ( [4.11| ). 
Superimposed to this curve the normalized variance computed numerically from the zeros 
of the Riemann zeta function is also shown. The agreement between the numerical results 
and the theoretical prediction is remarkable and the two curves are almost indistinguishable. 
Also displayed is the normalized GUE energy variance given by Eq. ( |3.12D . The saturation 
effect with respect to random matrix theory is clearly visible for a number of levels larger 
than n ~ 5. For higher values of n some oscillations are observed. These oscillations are 
well reproduced by the prime-number contributions ( [4.11|) . As predicted by the theory a 
damping of the oscillations is observable and illustrated in the inset in a wider range and 
on an expanded scale. The theory and numerical data are again in perfect agreement. For 
large values of n the curve tends to the constant 4.60, which is in very good agreement with 
the theoretical value given by Eq.(f4.9|), 1 — log Tmin = 4.56. 

It is also interesting to display as a function of n the behavior of the canonical energy 
variance for the zeros of the Riemann zeta function. As was mentioned in the intro- 
duction, this variance is not a two-point measure and therefore the semiclassical theories 
developed for the form factor do not apply. We could again expect that, even in the non- 
universal regime, for a large number of particles the canonical and grand canonical variances 
are very closely related. Fig. ^ shows that indeed this is the case. The normalized energy 
variance A^/(n^/27r^) presents the same non-universal oscillations in the saturation regime 
as Ve (shifted however to a lower average value). The non- universal features in imply via 
Eq.( |2.2| ) their presence in the autocorrelation function /(j). This is confirmed by a study of 
Odlyzko who has observed significant oscillatory deviations with respect to random matrix 
theory for large values of j 0, attributed to effects due to primes. 

As can be seen in Fig. ^ the behavior of the two variances A^ and Ve is quite different 
for n small. This difference is an (uninteresting) manifestation of the discreteness of n in 
the canonical case. 



15 



V. CONCLUDING REMARKS 



Assuming that the single particle spectrum is given by the eigenvalues of a random 
matrix, we have determined the typical fluctuations of the total energy of a system of non- 
interacting fermions. On the light of the original applications of random matrix theory 
this approach may look at flrst sight very unnatural. Indeed, these theories were originally 
applied to (nuclear) many body systems in spectral regions and for properties unrelated to 
the mean fleld, like for example neutron resonances of the compound nucleus. Here we are 
somehow adopting the opposite point of view, namely to start with an independent particle 
motion modeled by random matrix theory and neglect completely residual interactions. 
Surprisingly, there are physical situations for which this extreme view seems also to be 
fruitful, as should become clear by the end of this section. 

Two slightly different approaches have been investigated. In the flrst one, denoted 
"canonical", we consider the fluctuations of a flxed number n of occupied levels. These 
fluctuations are shown to be directly related to the autocorrelation function of consecutive 
spacings, for which we have proposed an ansatz. The parameters entering the ansatz have 
been determined from the asymptotic behavior of the fluctuations of the spacing variance 
a^(n). The canonical variance of the total energy is then computed, including corrections 
up to order one. The leading-order term is found to be log?T,/(/57r^). 

In the second - "grand canonical" - approach, the total energy variance of the single 
particle levels contained in an interval of given length has been considered. Contrary to the 
previous case, here the number of occupied levels is not flxed. We have computed the grand 
canonical variance in the limit of a large number of levels contained in the interval. Also an 
exact result has been given for /3 = 2. Both variances have the same leading-order behavior, 
but the higher order corrections are different. 

This difference in the higher order terms, and more generally the connection between 
"canonical" and "grand canonical" quantities, have been studied in detail and play a central 
role in the present study. Among the former we have considered the spacing variance cx^ and 
the energy variance A^, while the number variance T? and the energy variance Ve belong 
to the latter. The connection (|2.14| ) between the spacing and number variances has been 



essential in implementing the ansatz ( |2.7| ). As already discussed in 0, Eq.( |2.14|) , though 



asymptotic, remains a good approximation for any value of n. This is also conflrmed by our 
results. In fact, computing the difference to higher orders we flnd 

T?{n)-a^{n) = 1/6 + 0(1/^2) , 

with the coefficient of the 0{l/n^) term being a small (of the order 0.01) constant for any 
13. In contrast, for the other two quantities Eq. ( |3.10| ) is not a good approximation for small 
n. Computing the difference to higher orders we flnd 

V.W-A» 1 .;^_l(c_l:)\_l(,,_s,^^2^^o(llrf), (5.1) 



IT? 12 n 2 V 6/ n 3 r? 

where the constants have been determined in Section II. A. 

These relations between canonical and grand canonical quantities are those predicted 
by random matrix theory, and thus applicable in principle only in the universal regime (for 
example, Eq.(^]l]) describes the difference between the dot-dashed and the long dashed curves 
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in Fig. But surprisingly we find that Eqs.( |2.1^) and ( |5.1| ) always hold, irrespectively of 
the regime considered. This "universality in the non-universal regime" is illustrated in Fig. |^ 
where both differences for the zeros of the Riemann zeta function are plotted (as we have 
mentioned before, a transition to a non-universal regime appears for n > 5), and presumably 
it can be traced back to the validity of general incompressibility conditions for a fluid pT|,p|. 
An immediate and interesting consequence of this observation is that the non-universal 
corrections for canonical quantities (for which no explicit theory is available) ought to be 
identical to those of grand canonical quantities. The consequences of this remark and in 
particular the presence of non-universal contributions in the behavior of the /(i)'s are under 
study. 

The simplicity of periodic orbit theory in the interpretation of physical phenomena has 
demonstrated its power in different branches of physics. A well-known example is the pre- 
diction and the experimental confirmation in metallic clusters ||2^ of supershell effects 



in the density of states which can be basically understood in terms of a beating produced 
by two short periodic orbits of a spherical cavity having similar lengths. For irregularly 
shaped clusters, for which our model may be a starting point, the possibility of having al- 
most degenerate short periodic orbits is not excluded. However it is unclear how robust this 
situation may be under small perturbations, like for example the addition of particles to the 
system. As already mentioned, similar and perhaps more promising possibilities exist in the 
physics of quantum dots. 

Other systems for which the present study may be relevant concern diffusive systems 
in mesoscopic physics, like the motion of electrons in a disordered piece of metal. Due to 
the presence of impurities, in this case the dynamics of the electrons is not ballistic but 
rather described by a random walk. For short times the non-universalities are not related 
to short periodic orbits but to the return probability of a diffusive motion. This can be 
incorporated in the short time behavior of the form factor ||2^. In particular it does not 



lead to the saturation effect discussed here and presumably gives rise, for the variance of 
the total energy, to similar effects as the ones observed for the number variance which 
are governed by the diffusion constant and the space dimensionality. 

In nuclear physics there are some estimates of the fluctuations of the binding energies 
due to non-systematic effects. When the macroscopic part as well as shell corrections are 
substracted out from the measured total energy it is found that the r.m.s. of the remaining 
fluctuations is about 0.5 MeV for heavy nuclei Since the total binding energy for those 
nuclei is about 1000 MeV, the relative fluctuations due to non-systematic effects are of 
order R^xp = 0.5 x 10^^. With the present model, since A ^ n^/logn/n for /3 = 1 and 
< Xt >~ n'^/'2, we obtain a theoretical estimate of the relative fluctuations Rth ~ 7 x 10~^ 
for n ^ 200 which overestimates the experimental value by an order of magnitude. Before 
considering this disagreement as significant obvious effects neglected here like the energy 
variations of the mean level spacing should be taken into account. 
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APPENDIX. A 



Given a stationary sequence of successive random single particles energies . . . , Xj, Xj+i, . . ., 
take a point O at random on the real axis. It will lie between two levels, denoted by Xq and 
Xi- Call Si = Xi — Xq and d the (random) distance between xi and O. Fill the n successive 
levels located just above O and define the total energy Xl{n) with respect to O: 

n 

Xlin) = -i+l) Si +n d = X;{n - I) + n d , (A.l) 

i=2 

where Sj = Xi — Xi-i. By construction, d is the product of two independent random variables, 
e and si, with e distributed uniformly in the interval [0, 1] and si with the nearest neighbor 
spacing distribution of the original sequence. Using (esi) = (e) (si), (e) = 1/2 and (si) = 1 
one has 

(X;H) = (X,(n-l)) + ^ = ^. (A.2) 

To compute the variance A*^(n) of X/(n), care must be taken of the correlations between 
Si and the nearest neighbor spacings above Si. Using Eq. (|2.8| ) one obtains 



n n — 1 

(si Xt{n - 1)) = ^(n - z + 1) (si s,) = ^(^ " 3) Hj) = [c^Hn) - n ([s") - n)] /2 . 

(A.3) 

Following the same steps as in section II, using the fact that e is uncorrelated with si and 
Xt, using (e^sf) = {e^){s\) = (s^)/3, and using Eq. ( |A.3| ) one finally obtains 

A-(n) = A\n - 1) + y " ^) + f ^'H • (A.4) 



APPENDIX. B 

The energy variance is not independent of the reference point used to measure the energy. 
In the grand canonical case, when an arbitrary point a is used as reference point, f{x) in 
Eq.( |3.1|) should be taken as 

X — a —L < X < L /-n T\ 

|a:|>L. (^-^^ 

It follows from the definition of the energy variance ( |3.2| ) and ( |B.1| ) that the variance Veq 
with the reference point at the center of the energy interval considered and the energy 
variance Ve^ with the origin at related through the equation 

VeM = VE,in) + a^T?{L = n) , (B.2) 
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where T? is the number variance. This is an exact relation vahd for any statistical sequence 
of energy levels (uncorrelated, GE, etc). 

When computing the energy variance we have used as reference point for the energy the 
lower (or upper) end of the energy interval considered. Setting a = —L in Eq.( p.2| ), using 



Eq.( |3.9| ) and the asymptotic form of the number variance ( p.l2| ) we obtain for the variance 
with reference at the center of the interval 

VEoin) ^ -n'^logn + -i^- -rijn'^ , (B.3) 



with the constants given in Section II. A. The leading order term in Eq. (|B.3|) is twice smaller 



than the variance ( p.9| ) with the reference point at the lower (or upper) end of the interval. 

Because Veq and are positive definite, Eq.( p.2|) shows, incidentally, that the energy 
variance reaches a minimum when the reference point is at the center of the interval. 
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FIG. 1. Distribution function of the total energy Xt for n = 300. Histogram: from an ensemble 
of 500 Gaussian orthogonal matrices of dimension = 600; solid line: Gaussian distribution with 
mean value (Xt(300)) and variance A^(300) obtained from Eqs.(p]) and ( p.l7| ), respectively. 



22 




FIG. 2. Autocorrelation of spacings I{j) for (3 = 2. Dots: obtained from an ensemble of 500 
Gaussian unitary matrices of dimension = 500; solid line: the ansatz ( p. 71 ) for [5 = 2] squares: 
computed from the zeros of the Riemann zeta function (taken from [9]). 
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FIG. 3. Canonical number variance for /3 = 1 as a function of n. Solid line: from an 
ensemble of 500 Gaussian orthogonal matrices of dimension N = 600; dashed line: theory (see 
Eq.( 2.17D ); dot-dashed line: leading order term of the same equation. 
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FIG. 4. Normalized energy variances. Solid line: the normalized grand canonical variance 
VeI ir? I'i'T^'^^ for the zeros of the Riemann zeta function. In this curve the numerical results as 
well as the theoretical prediction are superimposed, and are almost indistinguishable. The same 
curve is plotted in a wider range and on an expanded scale in the inset. Dot-dashed line: random 
matrix result (/3 = 2) for the normalized grand canonical variance. Dashed line: the normalized 
canonical variance A^/ (n^/27r^) obtained numerically from the zeros of the Riemann zeta function. 
Long-dashed line: random matrix result (/3 = 2, obtained using the ansatz) for the normalized 
canonical variance. (See text for further explanation). 
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FIG. 5. Comparison between canonical and grand canonical quantities. The upper solid curve 
is the difference — cr^ obtained numerically from the zeros of the Riemann zeta function. The 
upper dashed horizontal line is at 1/6 ~ 0.166. The lower solid curve is the numerical result for 
the difference {Ve — A^)/n^ for the same set of zeros. The lower dashed curve displays Eq.( |5.1[ ), 
which tends to 1/12 ~ 0.0833 for large values of n. 
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